
#Simulation Replication

rm(list = ls()) #I know I shouldn't, but I'm doing it.

library(ggplot2)
library(gridExtra)
library(glue)
library(ggpubr)
library(scales)
library(fastDummies)
library(knitr)
library(ggpubr)
library(SimCorMultRes)
library(SimCorrMix)
library(simstudy)
library(stargazer)
library(survey)
library(NMOF)
library(devtools)
source_gist("https://gist.github.com/mrdwab/6424112")

#Setting up and running the simulations

countriesN <- 2000
uniqueID <- countriesN*1000


choices <- c(rep(1,210),rep(2,210),rep(3,160),rep(4,210),rep(5,210))


polProbs <- matrix(c(0.21, 0.21, 0.16, 0.21, 0.21,
                     0.21, 0.21, 0.16, 0.21, 0.21,
                     0.21, 0.21, 0.16, 0.21, 0.21),
                   nrow = 3, byrow = TRUE)

set.seed(25)

polDat <- genData(1000)
polOrdData <- genCorOrdCat(polDat, adjVar = NULL, baseprobs = polProbs, 
                   prefix = "q", rho = 0.5, corstr = "cs")

nonOrdData <- data.frame("id"=c(1:1000),"s1"=sample(choices,1000,replace =T), "s2"=sample(choices,1000,replace =T), 
                         "s3"=sample(choices,1000,replace =T))


combinedB <- merge(polOrdData,nonOrdData, by = "id")
combinedB$mq1 <- ifelse(combinedB$q1==3,1,0)
combinedB$mq2 <- ifelse(combinedB$q2==3,1,0)
combinedB$mq3 <- ifelse(combinedB$q3==3,1,0)


combinedB$ms1 <- ifelse(combinedB$s1==3,1,0)
combinedB$ms2 <- ifelse(combinedB$s2==3,1,0)
combinedB$ms3 <- ifelse(combinedB$s3==3,1,0)


testCensor <- (sum(combinedB$mq1)+sum(combinedB$mq2) +sum(combinedB$mq3))/3000 -
  (sum(combinedB$ms1)+sum(combinedB$ms2) +sum(combinedB$ms3))/3000

testCensor

combinedB$CostLying <- runif(1000, min = 0, max = 1)

CombinedFull <- combinedB[rep(seq_len(nrow(combinedB)), countriesN), ]
CombinedFull$country <- rep(1:countriesN, each=1000)

ExC1 <- c(runif(countriesN/2, min=0, max=0.1))
ExC2 <- c(runif(countriesN/2,min=0,max=1))

ExternalCost <- data.frame("country"=c(1:countriesN), 
                           "sensitivity"=c(ExC1,ExC2))

CombinedFull <- merge(CombinedFull,ExternalCost, by="country")

ExpressedQ11 <- c()
ExpressedQ2 <- c()
ExpressedQ3 <- c()

resulT11 <- -(CombinedFull$sensitivity)*(5-1)^2-CombinedFull$CostLying*(1 - CombinedFull$q1)^2
resulT12 <- -(CombinedFull$sensitivity)*(5-2)^2-CombinedFull$CostLying*(2 - CombinedFull$q1)^2
resulT13 <- -(CombinedFull$sensitivity)*(5-3)^2-CombinedFull$CostLying*(3 - CombinedFull$q1)^2
resulT14 <- -(CombinedFull$sensitivity)*(5-4)^2-CombinedFull$CostLying*(4 - CombinedFull$q1)^2
resulT15 <- -(CombinedFull$sensitivity)*(5-5)^2-CombinedFull$CostLying*(5 - CombinedFull$q1)^2
resultsJoin1 <- data.frame(resulT11,resulT12,resulT13,resulT14,resulT15)
resultsJoin1$Maxc <- max.col(resultsJoin1)
CombinedFull$rq1 <- resultsJoin1$Maxc

resulT21 <- -(CombinedFull$sensitivity)*(5-1)^2-CombinedFull$CostLying*abs(1 - CombinedFull$q2)^2
resulT22 <- -(CombinedFull$sensitivity)*(5-2)^2-CombinedFull$CostLying*abs(2 - CombinedFull$q2)^2
resulT23 <- -(CombinedFull$sensitivity)*(5-3)^2-CombinedFull$CostLying*abs(3 - CombinedFull$q2)^2
resulT24 <- -(CombinedFull$sensitivity)*(5-4)^2-CombinedFull$CostLying*abs(4 - CombinedFull$q2)^2
resulT25 <- -(CombinedFull$sensitivity)*(5-5)^2-CombinedFull$CostLying*abs(5 - CombinedFull$q2)^2
resultsJoin2 <- data.frame(resulT21,resulT22,resulT23,resulT24,resulT25)
resultsJoin2$Maxc <- max.col(resultsJoin2)
CombinedFull$rq2 <- resultsJoin2$Maxc

resulT31 <- -(CombinedFull$sensitivity)*(5-1)^2-CombinedFull$CostLying*abs(1 - CombinedFull$q3)^2
resulT32 <- -(CombinedFull$sensitivity)*(5-2)^2-CombinedFull$CostLying*abs(2 - CombinedFull$q3)^2
resulT33 <- -(CombinedFull$sensitivity)*(5-3)^2-CombinedFull$CostLying*abs(3 - CombinedFull$q3)^2
resulT34 <- -(CombinedFull$sensitivity)*(5-4)^2-CombinedFull$CostLying*abs(4 - CombinedFull$q3)^2
resulT35 <- -(CombinedFull$sensitivity)*(5-5)^2-CombinedFull$CostLying*abs(5 - CombinedFull$q3)^2
resultsJoin3 <- data.frame(resulT31,resulT32,resulT33,resulT34,resulT35)
resultsJoin3$Maxc <- max.col(resultsJoin3)
CombinedFull$rq3 <- resultsJoin3$Maxc

CombinedFull$rmq1 <- ifelse(CombinedFull$rq1==3,1,0)
CombinedFull$rmq2 <- ifelse(CombinedFull$rq2==3,1,0)
CombinedFull$rmq3 <- ifelse(CombinedFull$rq3==3,1,0)


CombinedFull$ComMiss <- CombinedFull$rmq1 + CombinedFull$rmq2 + CombinedFull$rmq3
CombinedFull$ComSensMiss <- CombinedFull$ms1 + CombinedFull$ms2 + CombinedFull$ms3

CombinedFull$AggFals <- abs(CombinedFull$rq1 - CombinedFull$q1) + abs(CombinedFull$rq2 - CombinedFull$q2) +
  abs(CombinedFull$rq3 - CombinedFull$q3)

testCensor2 <- (sum(CombinedFull$rmq1)+sum(CombinedFull$rmq2) +sum(CombinedFull$rmq3))/uniqueID -
  (sum(CombinedFull$ms1)+sum(CombinedFull$ms2) +sum(CombinedFull$ms3))/uniqueID
testCensor2

CountryData <- data.frame(aggregate(list("Political" = CombinedFull$ComMiss,
                                         "Sensitive" = CombinedFull$ComSensMiss, "AggFals"=CombinedFull$AggFals),
                                    by = list("country" = CombinedFull$country,"External"=CombinedFull$sensitivity), sum))
CountryData$IndexSens <- CountryData$Political/3000 - CountryData$Sensitive/3000
CountryData <- CountryData[order(CountryData$country),]
CountryData$Democracy <- c(rep("Democracy",countriesN/2),rep("Autocracy",countriesN/2))
CountryData$Bias <- CountryData$AggFals/3000
mean(CountryData$IndexSens)


#Figure 2
Fig1B <- ggplot(CountryData,aes(IndexSens, color = Democracy)) +
  geom_density() + xlab("Self-Censorship Index") + ylab("Density") +
  theme(legend.title = element_blank())

Fig1B
ggexport(Fig1B,filename = "Figure2.pdf")

#Figure 3
Fig2 <- ggplot(CountryData, aes(External,IndexSens, group = Democracy)) + 
  geom_smooth(aes(color=factor(Democracy)), se=F) + ylab("Self-Censorship Index") + xlab("External Costs") +
  theme(legend.title = element_blank()) 

Fig2a <- ggplot(CountryData, aes(Bias,IndexSens, group = Democracy)) + 
  geom_smooth(aes(color=factor(Democracy)), se=F) + ylab("Self-Censorship Index") + xlab("Sensitivity Bias") +
  theme(legend.title = element_blank()) 

Fig2b <- ggplot(CountryData, aes(External, Bias, group = Democracy)) + 
  geom_smooth(aes(color=factor(Democracy)), se=F) + ylab("Sensitivity Bias") + xlab("External Costs") +
  theme(legend.title = element_blank()) 

Fig2arr <- ggarrange(Fig2b,Fig2,Fig2a, common.legend = T, nrow = 1)
Fig2arr
ggexport(Fig2arr,filename = "Figure3.pdf")

#Figure 4

CountryExample <- CombinedFull[CombinedFull$country==27 | CombinedFull$country==1006 | CombinedFull$country==1012,]
CountryExampleSum <- data.frame(aggregate(list("Q1"=CountryExample$rq1, "Q2"=CountryExample$rq2, "Q3"=CountryExample$rq3),by = list("country" = CountryExample$country,"External"=CountryExample$sensitivity), sum))
CountryExample$obv <- 1
CountryExampleSum <- data.frame(aggregate(obv ~ country + sensitivity + rq1, CountryExample, sum),
                                aggregate(obv ~ country + sensitivity + rq2, CountryExample, sum),
                                aggregate(obv ~ country + sensitivity + rq3, CountryExample, sum))

CountryExampleSum$rq1 <- ordered(CountryExampleSum$rq1)
CountryExampleSum$rq1 <- ordered(CountryExampleSum$rq1,levels = c(5,4,2,1,3))
levels(CountryExampleSum$rq1)[5] <- "Don't know"
levels(CountryExampleSum$rq1)[1] <- "A great deal"
levels(CountryExampleSum$rq1)[2] <- "Quite a lot"
levels(CountryExampleSum$rq1)[3] <- "Not very much"
levels(CountryExampleSum$rq1)[4] <- "None at all"
CountryExampleSum$country <- factor(CountryExampleSum$country, levels = c(27,1012,1006))
levels(CountryExampleSum$country)[1] <- "Democracy"
levels(CountryExampleSum$country)[2] <- "Autocracy-1"
levels(CountryExampleSum$country)[3] <- "Autocracy-2"

CountryExampleSum$rq2 <- ordered(CountryExampleSum$rq2)
CountryExampleSum$rq2 <- ordered(CountryExampleSum$rq2,levels = c(5,4,2,1,3))
levels(CountryExampleSum$rq2)[5] <- "Don't know"
levels(CountryExampleSum$rq2)[1] <- "A great deal"
levels(CountryExampleSum$rq2)[2] <- "Quite a lot"
levels(CountryExampleSum$rq2)[3] <- "Not very much"
levels(CountryExampleSum$rq2)[4] <- "None at all"

CountryExampleSum$rq3 <- ordered(CountryExampleSum$rq3)
CountryExampleSum$rq3 <- ordered(CountryExampleSum$rq3,levels = c(5,4,2,1,3))
levels(CountryExampleSum$rq3)[5] <- "Don't know"
levels(CountryExampleSum$rq3)[1] <- "A great deal"
levels(CountryExampleSum$rq3)[2] <- "Quite a lot"
levels(CountryExampleSum$rq3)[3] <- "Not very much"
levels(CountryExampleSum$rq3)[4] <- "None at all"




FigCountry <- ggplot(CountryExampleSum, aes(y = (obv/10), x = country, fill = rq1)) + geom_bar(stat="identity") +
  xlab("") + ylab("Percentage") + theme(legend.title = element_blank()) +
  ggtitle("Question 1") 
FigCountry
FigCountry2 <- ggplot(CountryExampleSum, aes(y = (obv.1/10), x = country, fill = rq2)) + geom_bar(stat="identity") +
  xlab("") + ylab("Percentage") + theme(legend.title = element_blank()) +
  ggtitle("Question 2")
FigCountry2
FigCountry3 <- ggplot(CountryExampleSum, aes(y = (obv.2/10), x = country, fill = rq3)) + geom_bar(stat="identity") +
  xlab("") + ylab("Percentage") + theme(legend.title = element_blank()) +
  ggtitle("Question 3")
FigCountry3

tableCountry <- CountryExampleSum[c(5:7),c(CountryExampleSum$country, CountryExampleSum$sensitivity)]
tableCountry <- tableCountry[,c(1:2)]
tableCountry$SCI <- 0.042
tableCountry$Repression <- c("Low","Low","High")
tableCountry <- tableCountry[,c(1,2,4,3)]
FigCTable <- ggtexttable(format(tableCountry, digits=2), rows = NULL, 
                         cols = c("Regime","External Cost","Repression","SC Index"), 
                         theme = ttheme("minimal"))
FigCTable
FigCountryFull <- ggarrange(FigCountry,FigCountry2,FigCountry3, FigCTable,ncol = 2, nrow = 2)
FigCountryFull


ggexport(FigCountryFull,filename = "Figure4.pdf")


#MSCI
resulQ11 <- -(CombinedFull$sensitivity)*0.5*(4)^2 - CombinedFull$CostLying*(1 - CombinedFull$q1)^2
resulQ12 <- -(CombinedFull$sensitivity)*(0.5*(3)^2+0.5)-CombinedFull$CostLying*(2 - CombinedFull$q1)^2
resulQ13 <- -(CombinedFull$sensitivity)*(0.5*(2)^2+0.5*4)-CombinedFull$CostLying*(3 - CombinedFull$q1)^2
resulQ14 <- -(CombinedFull$sensitivity)*(0.5*(3)^2+0.5)-CombinedFull$CostLying*(4 - CombinedFull$q1)^2
resulQ15 <- -(CombinedFull$sensitivity)*0.5*(4)^2-CombinedFull$CostLying*(5 - CombinedFull$q1)^2
resultsJoinQ1 <- data.frame(resulQ11,resulQ12,resulQ13,resulQ14,resulQ15)
resultsJoinQ1$Maxc <- max.col(resultsJoinQ1)
CombinedFull$rqQ1 <- resultsJoinQ1$Maxc

resulQ21 <- -(CombinedFull$sensitivity)*0.5*(4)^2-CombinedFull$CostLying*abs(1 - CombinedFull$q2)^2
resulQ22 <- -(CombinedFull$sensitivity)*(0.5*(3)^2+0.5)-CombinedFull$CostLying*abs(2 - CombinedFull$q2)^2
resulQ23 <- -(CombinedFull$sensitivity)*(0.5*(2)^2+0.5*4)-CombinedFull$CostLying*abs(3 - CombinedFull$q2)^2
resulQ24 <- -(CombinedFull$sensitivity)*(0.5*(3)^2+0.5)-CombinedFull$CostLying*abs(4 - CombinedFull$q2)^2
resulQ25 <- -(CombinedFull$sensitivity)*0.5*(4)^2-CombinedFull$CostLying*abs(5 - CombinedFull$q2)^2
resultsJoinQ2 <- data.frame(resulQ21,resulQ22,resulQ23,resulQ24,resulQ25)
resultsJoinQ2$Maxc <- max.col(resultsJoinQ2)
CombinedFull$rqQ2 <- resultsJoinQ2$Maxc

resulQ31 <- -(CombinedFull$sensitivity)*0.5*(4)^2-CombinedFull$CostLying*abs(1 - CombinedFull$q3)^2
resulQ32 <- -(CombinedFull$sensitivity)*(0.5*(3)^2+0.5)-CombinedFull$CostLying*abs(2 - CombinedFull$q3)^2
resulQ33 <- -(CombinedFull$sensitivity)*(0.5*(2)^2+0.5*4)-CombinedFull$CostLying*abs(3 - CombinedFull$q3)^2
resulQ34 <- -(CombinedFull$sensitivity)*(0.5*(3)^2+0.5)-CombinedFull$CostLying*abs(4 - CombinedFull$q3)^2
resulQ35 <- -(CombinedFull$sensitivity)*0.5*(4)^2-CombinedFull$CostLying*abs(5 - CombinedFull$q3)^2
resultsJoinQ3 <- data.frame(resulQ31,resulQ32,resulQ33,resulQ34,resulQ35)
resultsJoinQ3$Maxc <- max.col(resultsJoinQ3)
CombinedFull$rqQ3 <- resultsJoinQ3$Maxc

CombinedFull$rmqQ1 <- ifelse(CombinedFull$rqQ1==3,1,0)
CombinedFull$rmqQ2 <- ifelse(CombinedFull$rqQ2==3,1,0)
CombinedFull$rmqQ3 <- ifelse(CombinedFull$rqQ3==3,1,0)


CombinedFull$ComMissQ <- CombinedFull$rmqQ1 + CombinedFull$rmqQ2 + CombinedFull$rmqQ3
CombinedFull$AggFalsQ <- abs(CombinedFull$rqQ1 - CombinedFull$q1) + abs(CombinedFull$rqQ2 - CombinedFull$q2) +
  abs(CombinedFull$rqQ3 - CombinedFull$q3)

CountryCounter <- data.frame(aggregate(list("Political" = CombinedFull$ComMissQ,
                                            "Sensitive" = CombinedFull$ComSensMiss, "AggFallsQ"=CombinedFull$AggFalsQ),
                                       by = list("country" = CombinedFull$country,"External"=CombinedFull$sensitivity), sum))
CountryCounter$IndexSens <- CountryCounter$Political/3000 - CountryCounter$Sensitive/3000
CountryCounter$Bias <- CountryCounter$AggFallsQ/3000
CountryCounter <- CountryCounter[order(CountryCounter$country),]
CountryCounter$Democracy <- c(rep("Democracy",countriesN/2),rep("Autocracy",countriesN/2))


#Figure A6

Fig3B <- ggplot(CountryCounter,aes(IndexSens, color = Democracy)) +
  geom_density() + xlab("Self-Censorship Index") + ylab("Density") +
  theme(legend.title = element_blank()) 

Fig3B
ggexport(Fig3B,filename = "FigureA6.pdf")


#Figure A4

Fig4 <- ggplot(CountryCounter, aes(External,IndexSens, group = Democracy)) + 
  geom_smooth(aes(color=factor(Democracy)), se=F) + ylab("MSCI Score") + xlab("External Costs") +
  theme(legend.title = element_blank()) 
Fig4

Fig4a <- ggplot(CountryCounter, aes(Bias,IndexSens, group = Democracy)) + 
  geom_smooth(aes(color=factor(Democracy)), se=F) + ylab("MSCI Score") + xlab("Sensitivity Bias") +
  theme(legend.title = element_blank()) 
Fig4a

Fig4b <- ggplot(CountryCounter, aes(External, Bias, group = Democracy)) + 
  geom_smooth(aes(color=factor(Democracy)), se=F) + ylab("Sensitivity Bias") + xlab("External Costs") +
  theme(legend.title = element_blank()) 
Fig4b

Fig4arr <- ggarrange(Fig4b,Fig4,Fig4a, common.legend = T, nrow = 1)
Fig4arr
ggexport(Fig4arr,filename = "FigureA4.pdf")


#Simulation Summary
combinedB$obv <- 1
SummarySimsDis <- data.frame(aggregate(obv ~ q1, combinedB,sum), aggregate(obv ~ q2, combinedB,sum), 
                             aggregate(obv ~ q3, combinedB,sum), aggregate(obv ~ s1, combinedB,sum),
                             aggregate(obv ~ s2, combinedB,sum), aggregate(obv ~ s3, combinedB,sum))


#Figure 1


indData <- CombinedFull[CombinedFull$id==14 |CombinedFull$id==60 | CombinedFull$id==814,]
indData <- data.frame("ID" = as.factor(indData$id),"True"=indData$q1,"Expressed"=indData$rq1,
                      "Internal"=indData$CostLying,"External"=indData$sensitivity)

indData$InternalF <- ordered(indData$Internal)
levels(indData$InternalF)[1] <- "Low"
levels(indData$InternalF)[2] <- "Medium"
levels(indData$InternalF)[3] <- "High"



FigureInd <- ggplot(indData, aes(x=External,y=Expressed, color=InternalF)) + geom_line(size=1.5) +
  xlab("Repression") + ylab("Expressed Preference") +
  guides(color=guide_legend(title="Internal Cost of Preference Falsification:")) +
  theme(legend.position = "bottom") + 
  scale_y_continuous(labels=c("Most Critical", "Critical","Don't Know", "Supportive","Most Supportive"))
FigureInd

ggexport(FigureInd,filename = "Figure1.pdf")


# T-Test differences (Footnote 7)
sdDemocracies <- sd(CountryData[CountryData$Democracy=="Democracy",]$IndexSens)
sdDemocracies
sdAut <- sd(CountryData[CountryData$Democracy=="Autocracy",]$IndexSens)
sdAut
testDifference <- t.test(CountryData$IndexSens~CountryData$Democracy, var.equal=F)
testDifference

#Table 1

tData <- data.frame(Response = c("None at all", "Not very much", "Don't know", "Quite a lot", "A great deal"),
                    q1 = as.vector(prop.table(table(combinedB$q1))), 
                    q2 = as.vector(prop.table(table(combinedB$q2))), 
                    q3 = as.vector(prop.table(table(combinedB$q3))),
                    q4 = as.vector(prop.table(table(combinedB$s1))), 
                    q5 = as.vector(prop.table(table(combinedB$s2))),
                    q6 = as.vector(prop.table(table(combinedB$s3))))
stargazer(tData,summary = F,
          covariate.labels = c("Response","Political Q1","Political Q2","Political Q3", "Non-Sensitive Q1",
                               "Non-Sensitive Q2", "Non-Sensitive Q3"), 
          out="Table1.tex")

#Table A1

dataForMatrix <- data.frame("RA1" = combinedB$q1, 
                        "RA2" = combinedB$q2, "RA3" = combinedB$q3,
                        "NS1" = combinedB$s1, "NS2" = combinedB$s2, "NS3" = combinedB$s3)
corMatrix <- cor(dataForMatrix)
stargazer(corMatrix, out = "TableA1.tex", digits = 2, title = "Corrrelation Matrix for Simulated Preferences")

#Figure A1

pdf("FigureA1.pdf")
plot(table(CountryData$IndexSens), ylab = "Frequency",xlab = "SCI Score")
dev.off()

#Figure A5
pdf("FigureA5.pdf")
plot(table(CountryCounter$IndexSens), ylab = "Frequency", xlab = "MSCI Score")
dev.off()

#Varying Support Levels
#Low Support for the Regime (High Opposition)

polProbsH <- matrix(c(0.31, 0.31, 0.16, 0.11, 0.11,
                      0.31, 0.31, 0.16, 0.11, 0.11,
                      0.31, 0.31, 0.16, 0.11, 0.11),
                   nrow = 3, byrow = TRUE)

polOrdDataH <- genCorOrdCat(polDat, adjVar = NULL, baseprobs = polProbsH, 
                           prefix = "q", rho = 0.5, corstr = "cs")

nonOrdDataH <- data.frame("id"=c(1:1000),"s1"=sample(choices,1000,replace =T), "s2"=sample(choices,1000,replace =T), 
                         "s3"=sample(choices,1000,replace =T))



combinedBH <- merge(polOrdDataH,nonOrdDataH, by = "id")
combinedBH$mq1 <- ifelse(combinedBH$q1==3,1,0)
combinedBH$mq2 <- ifelse(combinedBH$q2==3,1,0)
combinedBH$mq3 <- ifelse(combinedB$q3==3,1,0)

combinedBH$ms1 <- ifelse(combinedBH$s1==3,1,0)
combinedBH$ms2 <- ifelse(combinedBH$s2==3,1,0)
combinedBH$ms3 <- ifelse(combinedBH$s3==3,1,0)

combinedBH$CostLying <- runif(1000, min = 0, max = 1)

CombinedFullH <- combinedBH[rep(seq_len(nrow(combinedBH)), countriesN), ]
CombinedFullH$country <- rep(1:countriesN, each=1000)

CombinedFullH <- merge(CombinedFullH,ExternalCost, by="country")

resulT11 <- -(CombinedFullH$sensitivity)*(5-1)^2-CombinedFullH$CostLying*(1 - CombinedFullH$q1)^2
resulT12 <- -(CombinedFullH$sensitivity)*(5-2)^2-CombinedFullH$CostLying*(2 - CombinedFullH$q1)^2
resulT13 <- -(CombinedFullH$sensitivity)*(5-3)^2-CombinedFullH$CostLying*(3 - CombinedFullH$q1)^2
resulT14 <- -(CombinedFullH$sensitivity)*(5-4)^2-CombinedFullH$CostLying*(4 - CombinedFullH$q1)^2
resulT15 <- -(CombinedFullH$sensitivity)*(5-5)^2-CombinedFullH$CostLying*(5 - CombinedFullH$q1)^2
resultsJoin1 <- data.frame(resulT11,resulT12,resulT13,resulT14,resulT15)
resultsJoin1$Maxc <- max.col(resultsJoin1)
CombinedFullH$rq1 <- resultsJoin1$Maxc

resulT21 <- -(CombinedFullH$sensitivity)*(5-1)^2-CombinedFullH$CostLying*abs(1 - CombinedFullH$q2)^2
resulT22 <- -(CombinedFullH$sensitivity)*(5-2)^2-CombinedFullH$CostLying*abs(2 - CombinedFullH$q2)^2
resulT23 <- -(CombinedFullH$sensitivity)*(5-3)^2-CombinedFullH$CostLying*abs(3 - CombinedFullH$q2)^2
resulT24 <- -(CombinedFullH$sensitivity)*(5-4)^2-CombinedFullH$CostLying*abs(4 - CombinedFullH$q2)^2
resulT25 <- -(CombinedFullH$sensitivity)*(5-5)^2-CombinedFullH$CostLying*abs(5 - CombinedFullH$q2)^2
resultsJoin2 <- data.frame(resulT21,resulT22,resulT23,resulT24,resulT25)
resultsJoin2$Maxc <- max.col(resultsJoin2)
CombinedFullH$rq2 <- resultsJoin2$Maxc

resulT31 <- -(CombinedFullH$sensitivity)*(5-1)^2-CombinedFullH$CostLying*abs(1 - CombinedFullH$q3)^2
resulT32 <- -(CombinedFullH$sensitivity)*(5-2)^2-CombinedFullH$CostLying*abs(2 - CombinedFullH$q3)^2
resulT33 <- -(CombinedFullH$sensitivity)*(5-3)^2-CombinedFullH$CostLying*abs(3 - CombinedFullH$q3)^2
resulT34 <- -(CombinedFullH$sensitivity)*(5-4)^2-CombinedFullH$CostLying*abs(4 - CombinedFullH$q3)^2
resulT35 <- -(CombinedFullH$sensitivity)*(5-5)^2-CombinedFullH$CostLying*abs(5 - CombinedFullH$q3)^2
resultsJoin3 <- data.frame(resulT31,resulT32,resulT33,resulT34,resulT35)
resultsJoin3$Maxc <- max.col(resultsJoin3)
CombinedFullH$rq3 <- resultsJoin3$Maxc

CombinedFullH$rmq1 <- ifelse(CombinedFullH$rq1==3,1,0)
CombinedFullH$rmq2 <- ifelse(CombinedFullH$rq2==3,1,0)
CombinedFullH$rmq3 <- ifelse(CombinedFullH$rq3==3,1,0)

CombinedFullH$ComMiss <- CombinedFullH$rmq1 + CombinedFullH$rmq2 + CombinedFullH$rmq3
CombinedFullH$ComSensMiss <- CombinedFullH$ms1 + CombinedFullH$ms2 + CombinedFullH$ms3
CombinedFullH$AggFals <- abs(CombinedFullH$rq1 - CombinedFullH$q1) + abs(CombinedFullH$rq2 - CombinedFullH$q2) +
  abs(CombinedFullH$rq3 - CombinedFullH$q3)

CountryDataH <- data.frame(aggregate(list("Political" = CombinedFullH$ComMiss,
                                          "Sensitive" = CombinedFullH$ComSensMiss, "AggFals"=CombinedFullH$AggFals),
                                     by = list("country" = CombinedFullH$country,"External"=CombinedFullH$sensitivity), 
                                     sum))
CountryDataH$Bias <- CountryDataH$AggFals/3000
CountryDataH$IndexSens <- CountryDataH$Political/3000 - CountryDataH$Sensitive/3000
CountryDataH <- CountryDataH[order(CountryDataH$country),]
CountryDataH$Democracy <- c(rep("Democracy",countriesN/2),rep("Autocracy",countriesN/2))

#Figure A2

AFH2 <- ggplot(CountryDataH, aes(External,IndexSens, group = Democracy)) + 
  geom_smooth(aes(color=factor(Democracy)), se=F) + ylab("Self-Censorship Index") + xlab("External Costs") +
  theme(legend.title = element_blank()) 
AFH2

ggexport(AFH2,filename = "FigureA2.pdf")

#High Support for the Regime (Low Opposition)

polProbsL <- matrix(c(0.11, 0.11, 0.16, 0.31, 0.31,
                      0.11, 0.11, 0.16, 0.31, 0.31,
                      0.11, 0.11, 0.16, 0.31, 0.31),
                   nrow = 3, byrow = TRUE)


polOrdDataL <- genCorOrdCat(polDat, adjVar = NULL, baseprobs = polProbsL, 
                            prefix = "q", rho = 0.5, corstr = "cs")

nonOrdDataL <- data.frame("id"=c(1:1000),"s1"=sample(choices,1000,replace =T), "s2"=sample(choices,1000,replace =T), 
                          "s3"=sample(choices,1000,replace =T))



CombinedBL <- merge(polOrdDataL,nonOrdDataH, by = "id")
CombinedBL$mq1 <- ifelse(CombinedBL$q1==3,1,0)
CombinedBL$mq2 <- ifelse(CombinedBL$q2==3,1,0)
CombinedBL$mq3 <- ifelse(CombinedBL$q3==3,1,0)

CombinedBL$ms1 <- ifelse(CombinedBL$s1==3,1,0)
CombinedBL$ms2 <- ifelse(CombinedBL$s2==3,1,0)
CombinedBL$ms3 <- ifelse(CombinedBL$s3==3,1,0)

CombinedBL$CostLying <- runif(1000, min = 0, max = 1)

CombinedFullL <- CombinedBL[rep(seq_len(nrow(CombinedBL)), countriesN), ]
CombinedFullL$country <- rep(1:countriesN, each=1000)

CombinedFullL <- merge(CombinedFullL,ExternalCost, by="country")

resulT11 <- -(CombinedFullL$sensitivity)*(5-1)^2-CombinedFullL$CostLying*(1 - CombinedFullL$q1)^2
resulT12 <- -(CombinedFullL$sensitivity)*(5-2)^2-CombinedFullL$CostLying*(2 - CombinedFullL$q1)^2
resulT13 <- -(CombinedFullL$sensitivity)*(5-3)^2-CombinedFullL$CostLying*(3 - CombinedFullL$q1)^2
resulT14 <- -(CombinedFullL$sensitivity)*(5-4)^2-CombinedFullL$CostLying*(4 - CombinedFullL$q1)^2
resulT15 <- -(CombinedFullL$sensitivity)*(5-5)^2-CombinedFullL$CostLying*(5 - CombinedFullL$q1)^2
resultsJoin1 <- data.frame(resulT11,resulT12,resulT13,resulT14,resulT15)
resultsJoin1$Maxc <- max.col(resultsJoin1)
CombinedFullL$rq1 <- resultsJoin1$Maxc

resulT21 <- -(CombinedFullL$sensitivity)*(5-1)^2-CombinedFullL$CostLying*abs(1 - CombinedFullL$q2)^2
resulT22 <- -(CombinedFullL$sensitivity)*(5-2)^2-CombinedFullL$CostLying*abs(2 - CombinedFullL$q2)^2
resulT23 <- -(CombinedFullL$sensitivity)*(5-3)^2-CombinedFullL$CostLying*abs(3 - CombinedFullL$q2)^2
resulT24 <- -(CombinedFullL$sensitivity)*(5-4)^2-CombinedFullL$CostLying*abs(4 - CombinedFullL$q2)^2
resulT25 <- -(CombinedFullL$sensitivity)*(5-5)^2-CombinedFullL$CostLying*abs(5 - CombinedFullL$q2)^2
resultsJoin2 <- data.frame(resulT21,resulT22,resulT23,resulT24,resulT25)
resultsJoin2$Maxc <- max.col(resultsJoin2)
CombinedFullL$rq2 <- resultsJoin2$Maxc

resulT31 <- -(CombinedFullL$sensitivity)*(5-1)^2-CombinedFullL$CostLying*abs(1 - CombinedFullL$q3)^2
resulT32 <- -(CombinedFullL$sensitivity)*(5-2)^2-CombinedFullL$CostLying*abs(2 - CombinedFullL$q3)^2
resulT33 <- -(CombinedFullL$sensitivity)*(5-3)^2-CombinedFullL$CostLying*abs(3 - CombinedFullL$q3)^2
resulT34 <- -(CombinedFullL$sensitivity)*(5-4)^2-CombinedFullL$CostLying*abs(4 - CombinedFullL$q3)^2
resulT35 <- -(CombinedFullL$sensitivity)*(5-5)^2-CombinedFullL$CostLying*abs(5 - CombinedFullL$q3)^2
resultsJoin3 <- data.frame(resulT31,resulT32,resulT33,resulT34,resulT35)
resultsJoin3$Maxc <- max.col(resultsJoin3)
CombinedFullL$rq3 <- resultsJoin3$Maxc

CombinedFullL$rmq1 <- ifelse(CombinedFullL$rq1==3,1,0)
CombinedFullL$rmq2 <- ifelse(CombinedFullL$rq2==3,1,0)
CombinedFullL$rmq3 <- ifelse(CombinedFullL$rq3==3,1,0)

CombinedFullL$ComMiss <- CombinedFullL$rmq1 + CombinedFullL$rmq2 + CombinedFullL$rmq3
CombinedFullL$ComSensMiss <- CombinedFullL$ms1 + CombinedFullL$ms2 + CombinedFullL$ms3

CombinedFullL$AggFals <- abs(CombinedFullL$rq1 - CombinedFullL$q1) + abs(CombinedFullL$rq2 - CombinedFullL$q2) +
  abs(CombinedFullL$rq3 - CombinedFullL$q3)

CountryDataL <- data.frame(aggregate(list("Political" = CombinedFullL$ComMiss,
                                          "Sensitive" = CombinedFullL$ComSensMiss, "AggFals"=CombinedFullL$AggFals),
                                     by = list("country" = CombinedFullL$country,"External"=CombinedFullL$sensitivity), 
                                     sum))
CountryDataL$Bias <- CountryDataL$AggFals/3000

CountryDataL$IndexSens <- CountryDataL$Political/3000 - CountryDataL$Sensitive/3000
CountryDataL <- CountryDataL[order(CountryDataL$country),]
CountryDataL$Democracy <- c(rep("Democracy",countriesN/2),rep("Autocracy",countriesN/2))

#Figure A3

AFL2 <- ggplot(CountryDataL, aes(External,IndexSens, group = Democracy)) + 
  geom_smooth(aes(color=factor(Democracy)), se=F) + ylab("Self-Censorship Index") + xlab("External Costs") +
  theme(legend.title = element_blank()) 
AFL2
ggexport(AFL2,filename = "FigureA3.pdf")

#Figure 5

AFL3 <- ggplot() + geom_smooth(data = CountryDataL,aes(x=External,y=IndexSens,se=F,  color="red")) +
  geom_smooth(data = CountryData,aes(x=External,y=IndexSens,se=F, color = "blue")) +
  geom_smooth(data = CountryDataH,aes(x=External,y=IndexSens,se=F, color="black")) + 
  ylab("Self-Censorship Index") + 
  xlab("External Costs") +
  scale_color_manual(name = "True Support for the Regime", values =c("red"="red","black"="black","blue"="blue"), 
                     labels = c("Low", "Evenly Divided","High")) +
  theme(legend.position="bottom")
AFL3
ggexport(AFL3, filename = "Figure5.pdf")


